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Abstract 

We study the basic properties of accretion flows onto binary supcrmassivc black holes, including 
the cases in which a circumbinary disk is misaligned with the binary orbital plane, by means of three- 
dimensional Smoothed Particle Hydrodynamics simulations. We find that a circular binary system with a 
misaligned circumbinary disk normally produces a double peaked mass-accretion-rate variation per binary 
orbit. This is because each black hole passes across the circumbinary disk plane and captures gas twice 
in one orbital period. Even in misaligned systems, however, a single peaked mass-accretion-rate variation 
per binary orbit is produced, if the orbital eccentricity is moderately large (e > 0.3). The number of peaks 
in mass accretion rates can be understood simply in terms of the orbital phase dependence of the distance 
between each binary black hole and its closest inner edge of the circumbinary disk. In the cases of eccentric 
binary black holes having different masses, the less massive black hole can get closer to the circumbinary 
disk than the massive one, thus tidally splitting gas from its inner edge, but the created gas flows are 
comparably captured by both black holes with a short time delay. As a consequence, the combined light 
curve shows periodic occurrence of double-peaked flares with a short interval. This may account for the 
observed light variations of OJ287. 
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1. Introduction 

Hierarchical structure formation scenario tells us that a galaxy grows by merger of smaller galaxies. The recently 
confirmed correlation between the mass of supermassive black holes (SMBHs) and the mass or luminosities of the 
bulge of their host galaxies strongly support the idea that SMBHs have grown with the growth of their host galaxies 
(Magorrian et al. 1998; Gebhardt et al. 2000; Ferrarese & Merritt 2000). These relationships suggest that each SMBH 
in the center of each galaxy should have evolved toward coalescence in a merged galaxy. If this is the case, a binary of 
SMBHs should be formed in a merged galactic nucleus before two black holes finally coalesce, yet no binary SMBHs 
have clearly been identified so far except for some candidates at large separations (~ kpc) (e.g., Komossa et al. 2003). 

The merger of two SMBHs processes via three stages (Begelman et al. 1980): First, each black hole sinks into a 
common center of a merged galactic nucleus by the dynamical friction with the surrounding field stars and gas (Escala 
et al. 2005; Dotti et al. 2007). When the separation between the two black holes becomes as short as one parsec or 
so, the dynamical friction is no longer efficient and a hard binary is formed (Mayer et al. 2007). The separation of 
such a hard binary should have been reduced by some unknown mechanism. When the separation becomes as short 
as one parsec or less, finally, the binary rapidly merge by emitting gravitational wave radiation to become a single 
SMBH. However, there has been many discussions about by what mechanism the binary orbit decays in the second 
stage. One of predominant candidates is the interaction between the binary and the circumbinary disk (Ivanov et al. 
1999; Gould & Rix 2000; Armitagc & Natarajan 2002; Hayasaki 2009; Haiman et al. 2009; Cuadra et al. 2009; Lodato 
et al. 2009; Hayasaki et al. 2010; Nixon et al. 2011a; Kocsis ct al. 2012a; Kocsis et al. 2012b). 

It is an observational challenge how to identify binary SMBHs on such a subparsec scale (see Komossa (2006) 
for a review). Several ways have been proposed: Periodic optical and radio outbursts (e.g., OJ 287) (Sillanpaa, et 
al. 1988; Valtonen et al. 2011), wiggled patterns of the radio jet, indicating precessional motions on a parsec scale 
(Yokosawa & Inoue 1985; Lovanov&Roland 2005), X-shaped morphology of radio lobes (Merritt&Ekers 2002), double- 
peaked broad emission lines in active galactic nuclei (AGNs) (Gaskell 1996) and distant quasars (Bogdanovic et al. 
2009; Dotti et al. 2009; Boroson & Lauer 2009; Montuori et al. 2011; Montuori et al. 2011; Eracleous et al. 2012), 
double compact cores with a flat radio spectrum (Rodriguez et al. 2006), orbital motion of the compact core with a 
periodic flux variation (Sudou et al. 2003; Iguchi & Sudou 2010), and so on. 
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Fig. 1. Configuration of our model. There are two angles (/3,y) which specify the orientation of the circumbinary disk plane with 
respect to the binary orbital plane (x-y plane). 

With successive discoveries of binary black hole candidates and ongoing constructions of advanced gravitational 
wave detectors such as eLZSVl(Amaro-Seoane et al. 2012), much attention has been recently paid to electromagnetic 
signatures from binary SMBH systems, such as afterglows (Milosavljevic & Phinney 2005; Tanaka & Menou 2010; 
Tanaka et al. 2012), precursors (Chang et al. 2009; Bode et al. 2010; Hayasaki 2011; Farris et al. 2011; Farris et 
al. 2012; Bode et al. 2012), periodic emissions (Artymowicz & Lubow 1996; Lehto & Valtonen 1996; Hayasaki et al. 
2007; Hayasaki et al. 2008; Bogdanovic et al. 2008; MacFadyen & Milosavljevic 2008; Roedig et al. 2011; Sesana et al. 
2012; Shi et al. 2012; Noble et al. 2012; Daniel et al. 2012), and a dual-jet structure (Palenzuela et al. 2010; Mosta 
et al. 2010), in the context of a massive black hole coalescence. There are also several theoretical studies on the 
observational methodology to probe the presence of a gaseous disk around coalescing binary black holes through the 
waveform analyses of the emitted gravitational waves (Kocsis et al. 2011; Yunes et al. 2011; Hayasaki et al. 2012). 

In most of the previous studies, it has been assumed that the circumbinary disk is aligned with the binary orbital 
plane. However, the angular momentum vector of the binary does not always coincide with that of the circumbinary 
disk because the circumbinary disk would be formed independently of the formation of binary SMBHs. Therefore, 
the orientation of a circumbinary disk plane can be taken arbitrarily with respect to the binary orbital plane. It 
has been shown quite recently that the circumbinary disk which rotates in such an opposite direction as the binary 
rotation can be stable, if the misalignment angle between the binary orbital plane and the initial circumbinary disk 
plane is more than tt/2 (Nixon 2012), whereas the circumbinary disk which rotates in a prograde direction for the 
binary is stable, if the misalignment angle is less than tt/2. However, accretion flows onto binary SMBHs around the 
misaligned circumbinary disk rotating in a prograde direction have been poorly investigated. The misalignment angle 
is expected to produce complex light variations which never appear otherwise. Especially, two periodic outbursts per 
orbital period may be obtained, since there are two epochs in one orbital period, when the distance between the binary 
and the inner edge of circumbinary disk is shortest. By contrast, the binary orbital eccentricity produces a single, 
steep peak per binary orbit in a coplanar system (Hayasaki et al. 2007; Roedig et al. 2011). 

In this paper, we focus our study on basic properties of accretion flows onto binary SMBHs including the cases of 
a misaligned circumbinary disk. In the next section, we will describe our models and methods of calculations. The 
numerical results will be presented in section 3. In section 4, we provide a simple model in order to understand the 
numerical results. The final section is devoted to summary and discussion. 

2. Our models and calculation methods 

In this section, we first explain our models and next describe our calculation methods. 
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2.1. Initial settings 

Figure 1 illustrates a schematic picture of the setting of our model; binary black holes rotating each other are 
surrounded by a misaligned circumbinary disk. The binary is put on the x-y plane with its center of mass being at 
the origin. The masses of the primary and secondary black holes are represented by Mi and M2, respectively. We 
put a circumbinary disk around the origin. The unit vector of specific angular momentum of the circumbinary disk is 
expressed by (e.g. Pringle 1996) 

j A = (cos 7 sin /3, sin 7 sin /3, cos /3), (1) 

where (3 is the tilt angle between the circumbinary disk plane and the binary orbital plane, and 7 is the azimuth of 
tilt. The position vector of each black hole is given by 

n = (riCos0,rjsin</!>,O) [i = 1,2], (2) 

where 

(3) 

1 + ecoscp 

and 4> is the true anomaly (hereafter, we regard the true anomaly as the orbital phase) (e.g., Murray & Dermott 1999), 
771 = q/(l + q) and 772 = 1/(1 + q) with binary mass ratio q = M2/M1, and a = a\ + ai is the semi-major axis of the 
binary, where a\ = r\\a and 02 = r\ia. 

The circumbinary disk initially has a radially uniform density profile between the initial radius of the disk-inner 
edge, n n i, and the radius of r = 7"i n i + 0.05a (i.e., the initial width of the circumbinary disk is set to be 0.05a). The 
material in the circumbinary disk rotates around the origin on the circular orbit with the Keplerian rotation velocity. 
The vertical density structure of the circumbinary disk is exponential; i.e., we assume a hydrostatic balance with 
constant temperature in the vertical direction. The initial mass of the circumbinary disk is 1.0 x 10 -4 Mq. The disk 
temperature is assumed to be T = 30700K everywhere. Note that this corresponds to the typical central temperature 
of a standard disk at r = 2a around a single black hole with 10 8 M© for a given mass input rate of Mini = 1-0M© yr -1 
(Kato, Fukue & Mineshige 2008). 

In Tabic 1, wc summarize common model parameters for all models in the present study: total black hole mass 
(Mi +M2), semi-major axis of the binary orbit (a), initial mass of the circumbinary disk (A/disk), mass injection rate 
(Mi n j), disk temperature T^, and Shakura-Sunyaev viscosity parameter ass (Shakura & Sunyacv 1973). The relation 
between ass and the SPH artificial viscosity parameters will be described in section 2.3. The binary orbital period, 
Porb, is estimated to be approximately 9.4yr by using the Keplerian third law. 

Table 1. Common model parameters 



Total black hole mass Mi+M 2 =1O S M 

Semi-major axis a = 0.01 pc 

(Orbital period) (P orb = 27ra 3 / 2 /v/G(Mi + M 2 ) « 9.4 yr) 

Initial disk mass Mdisk = 1-0 X 10~ 4 M Q 

mass injection rate M m j = I.OMq yr -1 

Disk temperature T d = 30700 K 

Disk viscosity ass = 0.1 



2.2. Boundary conditions 

Gas particles are added to the outer edge of the circumbinary disk from its outside in an arbitrary angles at a 
constant rate of M- m j = I.OMq yr -1 . The inner edge of the circumbinary disk is determined by the balance between 
the tidal/resonant torque exerted by the binary black holes and the viscous torque of the circumbinary disk. We take 
an initial inner edge radius of the circumbinary disk to be rj n i = 2.5 for e = 0.5. In an equal- mass and circular binary, 
on the other hand, we take an initial inner edge radius as r m i = 1.68a, corresponding to the tidal truncation radius 
where the tidal torque of the binary equals to the viscous torque of the circumbinary disk (Papaloizou & Pringle 1977: 
see also Table 1 of Artymowicz & Lubow 1994). 

We set the outer calculation boundary at r = 6.0a, which is sufficiently far from the disk region so that the outer 
boundary should not affect the flow dynamics in the binary SMBH system. The SPH particles passing outward across 
the outer calculation boundaries are removed from the simulation box. 

The accretion radius depends on the mass of each black hole. The black holes are modeled by sink particles with 
the fixed accretion radius of r acc = 0.1a or 0.05a, depending on the black hole mass. Note that each accretion radius 
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is two orders of magnitude larger than the Schwarzshild radius of each black hole. Numerically, we remove all the 
particles that enter the region inside the accretion radius. 

2.3. Numerical method 

The simulations presented below are performed with the three-dimensional (3D) SPH code, which is based on a 
version originally developed by Benz 1990; Bcnz ct al. 1990; Bate ct al. 1995 and has been extensively used for various 
systems by many authors (e.g., Okazaki et al. 2002; Hayasaki et al. 2007; Okazaki et al. 2010; Takata et al. 2012). 
The SPH equations are composed of a mass conservation equation, a momentum equation with the SPH standard 
artificial viscosity, and an isothermal equation of state in substitution for an energy equation. These equations with 
the standard cubic-spline kernel are integrated using a second-order Runge-Kutta-Fehlberg integrator with individual 
time steps for each particle and a variable smoothing length (Bate et al. 1995), which results in saving enormous 
computational time when a large range of dynamical timescales are involved. 

The artificial viscosity commonly used in SPH consists of two terms: a term that is linear in the velocity differences 
between particles, which produces a shear and bulk viscosity, and a term that is quadratic in the velocity differences, 
which is needed to eliminate particle interpenetration in high Mach number shocks. The parameters o:sph and /3sph 
control the linear and quadratic terms, respectively. In the simulations shown in this paper, the artificial viscosity 
is adjusted so as to keep the Shakura-Sunyaev viscosity parameter agg =0.1 (Shakura & Sunyaev 1973), using the 
approximate relation ass = O.lasPH V# an d /%ph = (see Section 2.2 of Hayasaki et al. 2007), where h and H are 
the smoothing length of individual particles and the scale-height of the circumbinary disk, respectively. 

Sink particles (i.e. black holes) are orbiting around each other, following the Kcplerian third law, because pertur- 
bations by the interaction with SPH particles are negligible. In all models described in next subsection, total run 
time is 60 in a unit of P rb, and the simulation time t is normalized by P rb through this paper. The orbital phase is 
forwardly shifted by 0.03 at the end of run (t = 60) to correct the accumulation of numerical errors. 



Table 2. Individual model parameters. 



Model 


q = Ah/Ah 


e 


Tin/ a 


08,7) 


Al 


1.0 


0.0 


1.68 


(0,0) 


A2 


1.0 


0.0 


1.68 


(77/6,77/2) 


A3 


1.0 


0.0 


1.68 


(7r/4,7r/2) 


A4 


1.0 


0.0 


1.68 


(tt/6,0) 


Bl 


1.0 


0.5 


2.50 


(0,0) 


B2 


1.0 


0.5 


2.50 


(7r/6,7r/2) 


B3 


1.0 


0.5 


2.50 


(7r/4,7r/2) 


B4 


1.0 


0.5 


2.50 


(tt/6,0) 


CI 


0.5 


0.5 


2.50 


(0,0) 


C2 


0.5 


0.5 


2.50 


(7r/6,7r/2) 


C3 


0.5 


0.5 


2.50 


(7T/4,7T/2) 



2.4- Calculated models 

In the present study we are concerned with the observable quantities for various configurations of binary SMBHs 
with the circumbinary disk. Accordingly, we calculated eleven models, in total, by varying the binary mass ratio, 
binary orbital eccentricity, title angle, and azimuth of tilt. In Table 2, we summarize the model parameters: from 
the left to right, model number (first column), mass ratio of the secondary black hole to the primary one (second 
column), orbital eccentricity (third column), initial radius of the inner edge of the circumbinary disk in units of a 
(fourth column), and tilt angle and azimuth of tilt (/3,7) (fifth column, see Figure 1 and equation (1) for definitions 
of (3 and 7). 

While Model A is an equal-mass and circular binary, Models B and C are eccentric binaries with equal black hole 
masses (Models B) and unequal black hole masses (Model C). Through this paper, we assign Model A2 as a fiducial 
model. 

3. Accretion flows from a misaligned circumbinary disk to binary black holes 

In this section, we examine how the basic properties of accretion flows onto each black hole depend on /?, 7, e, and 
q, by performing 3D SPH simulations. 
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Fig. 2. Orbital phase dependence of mass accretion rates for all models. Models Al— A4, Models B1-B4, and Models CI— C3 are 
shown in the left panel, the middle panel, and the right panel, respectively. The data is folded on the orbital period over 40 < t < 60. 
Here, we redefine the orbital phase in order that the phase corresponds to the epoch of the periastron passage (see section 2.3). 
The binary black holes are at the periastron (or at the apastron) at phase 0.0 (0.5) for eccentric binaries (Models B1-C3). For 
clarity, we vertically offset the mass accretion rates for Model Al and Bl by +0.3, for Models A2, B2, and CI by +0.2, and for 
Models A3, B3, and C2 by +0.1 with respect to those of Models A4, B4, and C3, respectively. In the right panel, the solid line 
and dashed (red) line represent the mass accretion rates onto the primary black hole and that of the secondary one, respectively. 



Table 3. Summary of averaged mass accretion rates, averaged circularization radii, and the number of SPH particles at the end 
of the run (t=60) for all models. Since the mass is injected to the circumbinary disk at the rate of lMQ/yr - 1 , the averaged mass 
accretion rates indicate the fraction of the mass accretion rate onto each black hole to the mass injection rate. 



Model 



Al 

A2 
A3 
A4 



(Mi) 
[M yr- 



(M 2 ) 
[Moyi- 1 ] 



(r c ,x) 
[al 



<rc,2> 
[al 



0.045 
0.037 
0.028 
0.037 



0.045 
0.036 
0.028 
0.037 



0.033 
0.062 
0.055 
0.063 



0.033 
0.062 
0.056 
0.063 



SPH 



62157 
67774 
53866 
81149 



Bl 


0.018 


0.017 


0.038 


0.034 


61022 


B2 


0.012 


0.012 


0.059 


0.058 


90168 


B3 


0.007 


0.007 


0.051 


0.051 


55987 


B4 


0.007 


0.007 


0.053 


0.057 


58000 


CI 


0.013 


0.013 


0.037 


0.017 


74489 


C2 


0.014 


0.008 


0.056 


0.028 


54729 


C3 


0.007 


0.003 


0.056 


0.025 


51317 



3.1. Mass accretion rates onto binary black holes 

We first show in Figure 2 the orbital phase dependence of mass accretion rates for all models. Through this paper, 
the mass accretion rate is calculated by counting the number of SPH particles entering the accretion radius of each 
black hole per unit time. Each mass accretion rate is folded on the orbital period over 40 < t < 60. For eccentric 
binaries (Models B1-C3), the binary is at the periastron (or at the apastron) at phase 0.0 (0.5). For clarity, we 
vertically offset the mass accretion rates for Model Al and Bl by +0.3, for Models A2, B2, and CI by +0.2, and for 
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Models A3, B3, and C2 by +0.1 with respect to those of Models A4, B4, and C3, respectively. 

Hayasaki et al. (2007) already studied the aligned disk case (i.e., /? = 0) corresponding to Model Al. Since the mass 
continuously falls onto the binary from the circumbinary disk, there is no remarkable orbital-phase dependence in a 
circular orbit. As a result, no periodic variations are produced in neither mass accretion rates onto the two black 
holes nor the luminosity of accretion disks surrounding the two black holes. In other words, the accretion proceeds in 
a quasi-steady fashion, and, hence, no large light variations but small amplitude, random fluctuations are observed. 
More details are seen in Figures 2, 4 and 8 of Hayasaki et al. (2007). We should note, however, that even in the case of 
a circular binary, periodic variations may be produced by the formation of spiral patterns excited on the circumbinary 
disk (MacFadyen & Milosavljevic 2008; Daniel et al. 2012). Such features never appear in our calculations with limited 
spatial dimension of the circumbinary disk and, so, will be explored in future. 

The situations differ markedly, if the circumbinary disk is misaligned with respect to the binary orbital plane. In 
Models A2 and A3, the circumbinary disks are tilted by /? = ir/6 and 7r/4, respectively, while 7 is kept the same (i.e., 
7 = 7r/2). Both models produce mass-accretion-rate variations with (more than) two peaks in one orbital period, as 
shown in the left panel of Figure 2. This can be understood in terms of the geometrical effect. In these misaligned 
systems, two black holes stay away from the circumbinary disk plane in most of time, but twice per orbital period they 
cross the circumbinary disk plane. Thus, there are two chances in one orbital period that black holes approach the 
inner region of the circumbinary disk and strip gas from the circumbinary disk. This effect makes the mass accretion 
rates enhanced around the phases of 0.0 and 7r. In addition, there are two minor peaks in the mass accretion rate of 
Model A2. Their origin is, however, unclear at this moment. 

We also examine how the mass-accrction-ratc variations depend on the azimuth of tilt (7) in Model A4, in which 
we assign 7 = 0, keeping other parameters the same as those in Model A2. The resultant variation patterns of the 
mass accretion rate is the same as those in Model A2 but with a horizontal shift by 7r/4, as is expected. 

In Models B1-B4, more clear-cut periodic light variations in mass accretion rates are obtained as a natural conse- 
quence of the eccentric binary orbit with e = 0.5. The coplanar case (i.e., Model Bl) was also studied by Hayasaki et al. 
(2007), who found that binary SMBHs periodically emit intense pulses as a result of the periodic interaction between 
the binary and the circumbinary disk. The pulse period is therefore the orbital period. The variation amplitude is 
estimated to be a factor of 6 — 7 for Models Bl and B2, while it is a factor of 3 — 4 for Model B3 and B4. Obviously, 
interactions are strongest when the black holes are at the apastron, where the black holes are closest to the inner edge 
of the circumbinary disk. However, the mass accretion rate reaches a peak not at the phase (f> ~ 7r (near apastron 
passage) but at (f> ~ 0.0 (near periastron passage). This is because it takes about half an orbital period for the gas 
stripped from the circumbinary disk to freely fall onto the black holes. In Model B4, we assign 7 = 0, keeping other 
parameters the same as those in Model B2. The resultant mass accretion rate has a single peak per binary orbit, and 
the variation amplitude is smaller than that of Model B2. It is interesting to note that double-peak nature, which 
was observed in Models A2 and A3, is no longer noticeable in Models B2-B4. This is because the effect of orbital 
eccentricity is stronger than both the misalignment effect and the rotation effect in the azimuth direction, thus erasing 
the double-peak nature (see section 4). Note that it seems that the mass-accretion-rate variations have a double 
peaked structure at 4> ~ hi Model B2 and B3, but such a structure may be numerical artifacts. 

We finally performed simulations of Models C1-C3 to see what variations in mass accretion rates are produced by 
the binary with different masses. The results are shown for Models C1-C3 in the right panel, where the solid line 
and dashed (red) line represent the mass accretion rates onto the primary and secondary black holes, respectively. 
Single peaked light variations arc obtained for these models. The mass accretion rate onto the secondary black hole 
is slightly smaller than those onto the primary black hole. Remarkably, the peak phases slightly differ among the two 
black holes. This is understood as follows. The gas in the inner edge of the circumbinary disk is tidally stripped at the 
apastron when the secondary black hole gets closest to the circumbinary disk. The gas freely infalls onto the binary 
while sneaking around the secondary black hole. While a part of the gas accretes at around the periastron onto the 
primary black hole moving around the center of mass of the binary system, another part of the gas accretes onto the 
secondary black hole at the periastron passage. The resultant superposed light curve exhibit one double-peaked flare 
per binary orbit, and its interval between peaks is shorter than those in Models A2, A3, and A4. 

3.2. Accretion flow patterns 

In order to visualize how gas accretion onto binary SMBHs occurs, we show in Figure 3 sequences of surface density 
contours for Model A2 (upper panels) and Model B2 (lower panels) at four different phases; at the phases of t = 59.0, 
59.25, 59.5, and 59.75 from the left to the right, respectively. These panels are illustrated in the inertial frame, and 
both of the black holes and the circumbinary disk are rotating in the counterclockwise direction. Note that the black 
holes are rotating more rapidly than the gas in the circumbinary disk according to the Kepler's third law. The two 
small points represent the locations of the primary black hole and secondary black hole, respectively. The solid small 
circles surrounding them represent their accretion radii, which are set to be r acc = 0.1a from the center of each black 
hole. The dashed circles in both models represent the inner edge of circumbinary disk ~ 1.68a for Model A2 and the 
1:3 outer Lindblad resonance radius ~ 2.1a for Model B2 (see Artymowicz & Lubow 1994), respectively. These density 
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Fig. 3. A sequence of snapshots of the accretion flow from a circumbinary disk onto binary SMBHs in Model A2 (upper four 
panels) and Model B2 (lower four panels). The color contours of surface density are displayed in the chronological order from the 
left (a) to the right (d) for each model. The surface density is calculated by integrating density in the direction of the circumbinary 
disk plane. The two small points and solid circles denote the positions of two black holes and the accretion radii surrounding them, 
respectively. The figure is shown in the incrtial frame; i.e., both of the binary and the circumbinary disk are rotating around the 
center in the anti-clockwise direction. Note that the disk plane is tilted by /3 = n/6 from the binary orbital plane with 7 = it/ 2 in 
both models. The density levels of each panel can be seen in the color chart (—3.5 < logS < 1.5) at the right side of lower panel 
(d). The dashed circle represents the initial inner-edge radius r; n = 1.68a for Molde A2, while it is the 1:3 outer Lindbrad resonance 
radius ~ 2.1a for Moldcl B2. Annotated in each panel are the major scale in units of O.Olpc, time in units of P or b and the number 
of SPH particles -/VgPHi respectively. 



maps arc projected to the circumbinary disk plane (denoted as the x*-y* plane). 

First, we decribe the case of Model A2 (upper panels). Although the flow patterns shown in the upper four panels 
look quite similar to those in Model Al reported by Hayasaki et al. (2007), we notice some differences between them. 
At the elapsed time of t = 59.0 and 59.5 as shown in panels (a) and (c), the innermost part of the circumbinary disk 
is most strongly distorted by the tidal force of each black hole, and therefore the surface density of the innermost 
part of the circumbinary disk is enhanced at the two positions of ~ 37r/4 and ~ 77r/4 from the x* axis, where cusp 
structures are observed. This is because each black hole is located within the circumbinary disk plane at those times. 
By contrast, at the elapsed time of t = 59.25 and 59.75 as shown in panels (b) and (d), the cusp density structure 
disappears, though density enhancements are seen to some extent. This is because each black hole is most distant 
from the circumbinary disk plane. Now we understand that the time variations of the overall density distribution 
repeat every half orbital period, which is the period of the passage of the black holes across the circumbinary disk 
plane. This is responsible for the double-peaked variations in the mass accretion rates. 

Let us, next, examine the case of Model B2, which is illustrated in the lower four panels in Figure 3. At the phase of 
t = 59.25 as shown in panel (b), the separation between the two black holes is going to increase and, hence, black holes 
are approaching the inner edge of the circumbinary disk. At t = 59.5 the binary separation reaches its maximum and 
the distance to the circumbinary disk is at minimum (see panel (c)). The gas in the innermost part of the circumbinary 
disk is pulled out by the black hole from t = 59.25 to 59.5, thereby a tidal tail being formed and extending inward from 
the innermost part of the circumbinary disk. This tidal tail continuously grows, and at the phase of t = 59.75 as shown 
in panel (d) we see a bridge connecting the two tidal tails extending from the opposite sides of the circumbinary disk. 
Gas is supplied to the black hole in a next moment, which is clearly seen in panel (a). In fact, it takes about a half of 
binary period for gas to fall onto the black holes from the inner edge of the circumbinary disk. This is the reason why 
mass accretion rate reaches its maximum value at around phase zero. It should be noted that such a gas dynamical 
behavior is very similar to that in Model Bl, which was reported previously (Hayasaki et al. 2007). In other words, 
the effects of the tilt angle (3 is not so appreciable for high eccentricity cases. We will discuss why it is so in section 4. 
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Fig. 4. Schematic view explaining the definitions of the two angles (4>,9) and the descending node. The orbital phase <f> and the 
azimuth angle 9 are measured from x-axis and the descending node, respectively. 

3.3. Averaged mass accretion rates and circularization radii 

In Table 3, we summarize the averaged mass-accretion rates and averaged circularization radii around black holes, 
and the number of SPH particle at the end of run for all models, where the circularization radius, r c ,i, is defined by 
r c,i = jt /GMi, [i — 1,2] and ji represents the specific angular momentum of the SPH particles which enter inside 
each accretion radius. Thus, their specific angular momentum is proportional to 1/2 power of circularization radius. 

It is noted from Table 3 that the same amount of mass injection is assumed for Models A1-A4, but the nevertheless 
the averaged accretion rate is higher in Model Al than others. These trends are also seen in Models B and C. 

The averaged circularization radii are about two orders of magnitude larger than the Schwarzshild radius corre- 
sponding to each black hole mass. This suggests that the averaged circularization radius of infalling material indicates 
the size of an accretion disk formed around each black hole. We also note that the averaged circularization radii for 
the case of j3 ^ is larger than those for the case of j3 = in all models. This is because the distance between each 
black hole and its nearest inner edge of the circumbinary disk is longer in the misaligned system (/3 ^ 0) than that of 
the coplanar system (/3 = 0). The mass tidally stripped from the circumbinary disk does not directly accrete onto each 
black hole but via the accretion disk around each black hole. We will briefly discuss how the accretion disk evolve in 
section 5. 

4. Simple semi-analytical models 

To understand in a simpler way how a variety of variation patterns in mass accretion rates arises, we construct a 
semi-analytical model. Since mass accretion occurs by tidal stripping of gas from the inner edge of the circumbinary 
disk, and since the gravitational attraction force to the circumbinary disk is strongest when the distance between the 
black hole and its nearest inner edge of the circumbinary disk is at a minimum, we can guess the number of peaks per 
orbital period by calculating the minimum distance as a function of the orbital phase for parameter sets of all models. 

The position vector of the inner edge of the circumbinary disk, r; n , can be expressed by 

fin = (ri n [cos6>sin7 + sin6'cos7COS/3],rj n [sin#sin7COs/3 — cosf?cos7], — ri n sin0sin/3), (4) 

where the orbital phase <fi and the azimuth angle 9 are measured from x-axis and the descending node, respectively 
(see Figure 4). From equations (2) and (4), we obtain the formula for the distance: 



d(/3,7,e,g;0) = |r in — = J 4a 2 + rf — 4ari[sin0cos(0 — 7)cos/3 — cos#sin(0 — 7)], 



(5) 
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where the value of 9 for a given (f> is numerically chosen so as to give the minimum distance. Here, we assign i = 2, 
since the less massive black hole can move on a larger extent, thereby getting closer to the circumbinary disk than the 
massive one. 

Figure 5 shows the orbital phase dependence of the distance between the black hole and its nearest inner edge of 
circumbinary disk. There are four different parameters: orbital eccentricity (e), mass ratio (q), tilt angle (/?), and 
azimuth of tilt (7). The dependence of d(/3,j,e,q;(f>) normalized by the semi- major axis on these parameters are shown 
in panels (i)-(v). The fiducial parameters are those of Model A2; i.e., (/3,7,e,g) = (7r/6,7r/2,0.0, 1.0). 

Panel (i) shows 7-dependence of normalized distance, d((3,j,e,q;(j>)/a, for fixed values of (q,e,(3) = (1.0,0.0,7r/6). The 
solid line, dotted line, dashed line, and dash-dotted line represent d((f>)/a for 7 = tt/2, 0, tt/6, and 7r/4, respectively. 
We note that d((f>)/a is the shortest twice per binary orbit and that the phases of the minimum distance shift by 
varying the parameter 7. Here, let us go back to see the difference between mass accretion rate of Model A2 and that 
of Model A4. The phase difference between the solid line (7 = tt/2) and the dotted line (7 = 0) is tt/2. This supports 
our simulation results of Models A2 and A4. 

Panel (ii) shows /3-dependence of the normalized distance for the fixed values of (q,e,j) = (1.0,0.0,7r/2). The solid 
line, dashed line, dash-dotted line, and dotted line represent d(<fi)/a for f3 = tt/6, 0, 7r/4, and tt/2, respectively. The 
normalized distance reaches its minimum value twice per binary orbital period except for the case with f3 = 0, in 
which d((f>)/a is constant. That is, the non-zero values of (3 is essential to produce two peaks per binary orbit in the 
mass-accretion-rate variations. This panel provides a reasonable explanation about why the results of Models A2-A4 
are shown in the left panel of Figure 2. 

Panel (iii) shows e-dependence of the normalized distance for the fixed values of (q,/3,7) = (1.0,7r/6,7r/2). The solid 
line, dashed line, dash-dotted line, dotted line, and three-dotted line represent d((j))/a for e = 0, 0.1, 0.2, 0.3, and 0.4, 
respectively. The normalized distance has two round peaks for e = 0, but the larger e is, the deeper becomes a hollow 
at around <fi = tt. Although we have mentioned that non-zero values of (3 give rise to the double peaked modulations 
in mass accretion rates, this is not always the case when the orbital eccentricity is not zero. This will explain that the 
mass accretion rates of Models B2, B3, and B4 have a single peaked shape. We have also calculated the cases with 
7 = 0.0, reaching the same conclusion; i.e., single peaked variations arc found for e > 0.3. This supports the result of 
Model B4. 

Panel (iv) shows g-dependence of the normalized distance between the secondary (less massive) black hole and its 
nearest inner edge of the circumbinary disk for the fixed values of (e,/3,7) = (0,tt/6,tt/2). The solid line, dashed line, 
dash-dotted line, dotted line, and three-dotted line represent d((j>)/a for q = 1.0, 0.1, 0.3, 0.5, and 0.7, respectively. The 
lower the binary mass ratio is, the shorter becomes d(<f)/a. Variations in the mass ratio do not change the number of 
peaks of the normalized distance. Obviously, the smaller q is, the smaller becomes the mean distance. This is because 
the semi-major axis of the secondary black hole (r-i in equation 3) increases as q decreases. 

Panel (v) shows the same 7-dcpendcncc as those in panel (i) but for e = 0.1. This panel shows how the normalized 
distance changes with 7 for the cases with an even smaller orbital eccentricity. The solid line, dotted line, dashed line, 
and dash-dotted line represent d((j>)/a for 7 = tt/2, 0, tt/6, and tt/4, respectively. The solid line is the same as the 
dashed line of panel (iii) . It is interesting to note that double peaked variation curve for the case of 7 = changes to 
a single peaked one as 7 increases. 

5. Discussion 

Periodic light variations are not generally observed from single SMBH systems. There is a report for the presence of 
quasi-periodic oscillations (QPOs) in AGNs (Gierlihski et al. 2008), but the QPOs do not mean periodic occurrence of 
clear-cut flares but mean rather gradual variations. Unlike single SMBH systems, periodic light variations can be seen 
in an eccentric binary black hole system because of periodically enhanced binary-disk interactions. In most of previous 
studies, it was implicitly assumed that the binary orbital plane is aligned with the circumbinary disk plane. Since 
there are no strong reasons to believe the binary-disk alignment, it is natural to relax this assumption and to examine 
what light variations are expected for the misaligned system. As shown in Figure 2, a circular binary surrounded by 
a misaligned circumbinary disk exhibits a double peak structure in the variations of mass accretion rates, whereas an 
eccentric binary with a moderately large orbital eccentricity shows a single peak per orbit, even if the circumbinary 
disk is inclined from the binary orbital plane. 

5.1. Inner edge radius of misaligned circumbinary disks 

The radius of the inner edge of the circumbinary disk is determined by the balance between the viscous torque of the 
prograde circumbinary disk and the tidal/resonant torque which acts on it (Artymowicz & Lubow 1994). In the case 
of moderate orbital eccentricity, the typical inner edge radius is at the 1 : 3 outer Lindblad resonance radius, which is 
estimated to be ~ 2.1a. If the circumbinary disk is misaligned with respect to the binary plane, the resonant torque 
will be weaker than otherwise, which makes the inner edge smaller. In order to more precisely determine the size of 
the inner edge of the circumbinary disk, we need to investigate how the tidal/resonant torque acts on the misaligned 




Fig. 5. Orbital phase dependence of the normalized distance, d(/3, 7, e,q; <f>) /a, between the black hole and the inner edge of 
circumbinary disk. The dependence of the normalized distance on the four parameters: 0, 7, e, and q arc shown in panels 
(i)— (v). In all the panels, the solid (black) line indicates the fiducial model; (0,y,e,q) = (tt/6, tt/2, 0.0, 1.0) for panels (i)— (iv) 
and (/3,7,e,q) = (rc/6, n/2, 0.1, 1.0) for panel (v). For the meanings of other lines, see text and the explanations in each panel. 



circumbinary disk. This is also essential to calculate the circumbinary disk structure. 

5.2. Triple disk model for 03 287 

A blazer OJ287 exhibits quasi-periodic optical outbursts of 12 year interval. In order to explain such quasi-periodic 
nature, a binary black hole model was firstly proposed by (Sillanpaa, et al. 1988). Lehto & Valtonen (1996) subsequently 
proposed the modified model that a less massive black hole orbits around a more massive one with undergoing 
relativistic precession and impacts twice per binary orbit on the accretion disk around the more massive black hole. 
This model can explain the quasi-periodicity of outbursts with a double peak structure at a few year interval. 

As shown in Models CI and C2 at the right panel of Figure 2, the solid line and dashed (red) line represent the 
mass accretion rates of the primary black hole and that of the secondary black hole, respectively. Respective mass 
accretion rates periodically vary with time and have a single peak around the periastron but slightly different orbital 
phases. The resultant combined mass accretion rate provides a periodic double peaked structure with a short interval. 
We, therefore, propose an alternative binary black hole model that the circumbinary disk around an eccentric binary 
composing of black holes with different masses is misaligned by a relatively small tilt angle (0 < /3 < tt/6). 

5.3. Precessions of misaligned circumbinary disks 

The binary-disk interaction also gives rise to a precession of the circumbinary disk (MacFadyen & Milosavljevic 
2008; Hayasaki & Okazaki 2009; Daniel et al. 2012). With the assumption that r»a and M2 <C Mi, the precession 
frequency of the misaligned circumbinary disk is given by (e.g., Nixon et al. 2011b) as 

Pprec Al + q ( r\V2 \ 



prcc 
orb 



3 q \aJ |cos/3| 



where P pre c and P or b are the precession period and binary orbital period, respectively. It is clear that P pr0 c is much 
longer than P rb, e.g. Pp roc ^ 100P O rb for q = 0.1 and r = 2a. This can produce a light variation in circumbinary disk 
with a beat period Pbeat: where it is expressed by 1/Pbeat = 1/Porb — 1 /Pprec- From equation (6), we obtain the beat 
period normalized by the binary orbital period: 
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which is slightly longer than the orbital period, as far as r > a. The beat period for the case of q = 0.1 and r = 2a is, 
for example, approximately 1.01P O rb- 111 a subsequent paper, we will study this topic in more detail. 

5.4- Luminosity variations 

Another interesting topic is how the mass finally accretes onto each black hole via an accretion disk from the 
circumbinary disk. In section 3.3, we show that the infalling material is circularized around each black hole (the 
estimated circularization radii can be seen in Table 3). This leads to the formation of the triple disk system, which 
is composed of two accretion disks around black holes and one circumbinary disk surrounding them (Hayasaki et 
al. 2008). However, it is poorly known what structure each accretion disk has and how it evolves. Once these two 
accretion disks are formed by gas supply from the circumbinary disk, they viscously evolve and gas in the accretion 
disk finally accretes onto each black hole after the viscous timescale. Assuming that the accretion disk is the standard 
disk, the viscous timescale measured at the circularization radius is much longer than the binary orbital period. The 
precise shape of light curves may be different from that of variations in mass accretion rates. It is interesting to 
examine the basic properties of radiations emitted from each accretion disk. 

6. Conclusions 

We have carried out numerical simulations of accretion flows from a circumbinary disk which is inclined from the 
binary orbital plane, in order to examine to what extent the basic properties of mass-accretion-rate variations may 
alter, compared with the coplanar cases. Our main conclusions are summarized as follows: 

1. We find that the mass accretion rates exhibit a double peak per binary orbit in a circular binary system, when 
the circumbinary disk is misaligned with the binary orbital plane. This is because each black hole passes across 
the circumbinary disk plane twice per binary orbit and then attracts the gas there. This double peak nature 
of mass accretion rates is also independent of the azimuth of tilt. The tilt angle is one of important orbital 
parameters to determine variation patterns of radiations emitted from a binary black hole system. 

2. The orbital eccentricity remains to be an important orbital parameter to produce single sharply peaked variations 
per binary orbit in mass accretion rates even in a misaligned circumbinary disk system. This is because each 
black hole is closest to the inner edge of circumbinary disk once per binary orbit in most cases. The simple 
semi-analytic model (see section 4) predicts that this single peak nature is independent of both the tile angle 
and the azimuth of tilt, as long as 0.3. 

3. In the case of an eccentric binary composing of black holes with different masses, the less massive black hole 
can get closer to the circumbinary disk than the massive one, thus tidally splitting gas from its inner edge, but 
the created gas flows are comparably captured by both black holes with a short time delay. The superposed 
accretion rates show periodic outbursts with an apparent double peaked structure with a short interval. 
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